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ABSTRACT 

We investigate supersonic, axisymmetric magnet ohydro dynamic (MHD) jets with a time-dependent 
injection velocity by numerical simulations with the PLUTO code. Using a comprehensive set of 
parameters, we explore different jet configurations in the attempt to construct models that can be 
directly compared to observational data of microjets. In particular, we focus our attention on the 
emitting properties of traveling knots and construct, at the same time, accurate line intensity ratios 
and surface brightness maps. Direct comparison of the resulting brightness and line intensity ratios 
distributions with observational data of microjets shows that a closer match can be obtained only 
when the jet material is pre- ionized to some degree. A very likely source for a pre- ionized medium is 
photoionization by X-ray flux coming from the central object. 

Subject headings: ISM: jets and outflows - (ISM): Herbig-Haro objects - Magnetohydrodynamics 
(MHD) - Shock waves - Methods: numerical 
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1. INTRODUCTION 

Jets from Young Stellar Objects (YSOs) derive their 
emission from the gas that has been heated and com- 
pressed by shocks. In fact, the actual jet matter is 
invisible for most of its extension and only the cool- 
ing zones behind the shocks emit a variety of lines 
that can be revealed with great accuracy and are 
rich of diagnostic indications on the post-shock phys- 
ical parameters such as temperature, density, ioniza- 
tion fraction and radial velocity. We refer especially 
to the so -called "microjets" like HH 30, DG Tau and 
RW Aur (IBacciotti et all 1 19991: lHa rtigan fc Morsel l2007t 
IBacciotti et all 12002: Me lnikov et al.l 12009( 1. where the 
line emission is limited to a region of the jet going up to 
about 4" — 5" from the forming star. Therefore, a careful 
study of the shock formation and evolution is crucial to 
understand the physical processes at work and to con- 
strain jet parameters that cannot be directly observed, 
such as the magnetic field intensity, the pre-shock density 
and temperature and the jet velocity. 

Radiative shocks have been studie d in steady-state 
cond iti ons by several authors (e.g., iCox fc Raymond! 
I1985L lHartigan et al.l Il994f h who derived the one- 
dimensional post-shock behavior of various physical pa- 
rameters (temperature, ionization fraction, electron den- 
sity, etc.) as functions o f the distance fr o m the shock 
front. More r ecently , Massagl ia et al.l (|200 5a) and 
Tesil eanu et al.l (|2009al ) have carried out ID numerical 
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studies of the time-dependent evolution of radiating, 
magnetized shocks. They have applied the results to the 
cases of DG Tau and HH 30, with the goal to reproduce 
the observed behavior of the line intensity ratios along 
the jet. 

These studies, as discussed by iRaga et all (|2007[ ) as 
well, brought about the problem of the numerical res- 
olution that is needed for a correct treatment of the 
post-shock region, especially as far as the reproduction 
of the line ratios is concerned. To solve this problem 
the authors have employed Adaptive Mesh Refinement 
(AMR) techniques, that allows to follow with great ac- 
curacy the sudden temperature drop behind the shock 
fr ont and save computati onal time. 

Tesilea nu et al.l (|2009af ) discussed as well the influence 
on the results of the cooling function details. They con- 
cluded that the use of a detailed treatment of radiative 
emissions and ionization/recombination processes in the 
jet material, as well as adequate numerical resolution are 
very important for the reproduction of emission line ra- 
tios, which are extremely sensitive parameters. Instead, 
to describe the general morphology of the jet and inte- 
grated emission line luminosities, an approximation of 
the total radiative losses gives good results, provided the 
numerical reso lution suffice to minimize numerical dis si- 
pation effects (|Raga et all 120071 : ITesileanu et alJl2Q0l . 

Even though these results were obtained in the ID 
limit, nonetheless they can serve as a guideline for multi- 
dimensional case, where additiona l physical effects, such 
as rotation feacciotti et al. 2002), can aff ect the shock 
evolution. ITesileanu et al.l (j2009b[ ) and iMignone et al.l 
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([20091 ) have carried out preliminary studies of the evo- 
lution of 2D shocks traveling along a jet deriving syn- 
thetic emission maps, full synthetic spectra and position- 
velocity (PV) diagrams of single lines. 

Previous theoretical investigations (e.g., 
Shull fc McKed (|1979h ) focused on stationary shock 
models where strong shocks were able to produce, 
via the UV radiation emission, the ionization of the 
pre-shock material, affecting the emission properties . 
Another approach was the one of (Hartiga n et al.l fl987). 
that provided for the interpretation of observational 
data a set of plane-parallel shock models, including some 
with totally ionized pre-shock medium, with the relative 
emission line fluxes. It was noticed, at that time, that 
large differences in the emission properties are related 
to the pre-ionization state of the pre-shock medium. 

In this paper, we consider the axisymmetric evolution 
of a train of shocks as they trave l alo ng; a jet, differ- 
ently fr om iMassaglia et al.l (|2005af ) and iTesileanu et al.l 
(2009a) that studied a single shock. 

These shocks are produced by imposing a sinusoidal 
perturbation on the jet structure, otherwise in radial 
equilibrium with the external environment. 

The use of AMR allows a careful treatment of the post- 
shock region, providing (at the highest level of refine- 
ment) a minimum grid size corresponding to about 0.02 
AU. The emissivity distribution o btained by numerical 
simul ations with the PLUTO code ([Mignone et al.l l2QQ7. 
l2011f ) is convolved with a point-spread-function (PSF) 
similar to the one of the observing instruments for com- 
parison with observations. As we shall see, a substantial 
improvement in reproducing the observed emission fea- 
tures can be achieved by introducing a pre-ionizati on of 
the jet material. Indeed, as recently pointed out ([Giidell 
regions surrounding proto-stars are subject to 
the action of X-rays able to ionize jet material to an im- 
portant degree that, due to the low recombination rates, 
lasts up to large distances from the jet origin. 

The plan of the paper is the following: In Section 
2 we discuss the initial equilibrium, perturbation, pre- 
ionization and parameters and the adopted techniques to 
model the problem; in Section 3 we present the results 
for different choice of the parameters; the conclusions are 
drawn in Section 4. 



2. THE MODEL 

Our model consists of a stationary jet model with a su- 
perimposed time-dependent injection velocity that pro- 
duces a chain of perturbations eventually steepening into 
shock waves. In what follows, the fluid density, velocity, 
magnetic field and thermal pressure will be denoted, re- 
spectively, with p, V (Vr.V^.Vz), B = (Br.B^.Bz) 
and p. The gas pressure depends on the plasma density 
p, temperature T and composition through the relation 
p = pksT/ (funn), where \i is the mean molecular weight 
and ks is the Boltzmann constant. 

Simulations are carried out by solving the time- 
dependent MHD equations in cylindrical axisymmetric 
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where p is the mass density, v = (y r , v<f>, v z ) is the veloc- 
ity, B = (B r , B^^Bz) the magnetic field, p t = p + B 2 /2 
denotes the total pressure, £ = — v x B is the electric 
field and E the total energy density: 
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with T = 5/3 the specific heat ratio. Also, = pViVj — 
BiBj are the flux dyad components and q_ = (E + p t )v — 
B(v • B) is the energy density flux. The source term Se 
accounts for radiative losses and is directly coupl e d to t he 
ionization network described in Tesile anu et al.l ((2008), 
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where k and i identify the element and its ionization 
stage, respectively, and S K j is a source term accounting 
for ionization and recombination processes. Given the 
range of temperature and density, we include the first 
three ionization stages of C, O, N, Ne, S besides hydrogen 
and helium. 

Numerical simulations have been performed in the 
computational domain defined by r £ [0,400] and z G 
[0, 1200] AU covered by a base grid of 128 x 384 cells, with 
6 additional levels of refinement with consecutive grid 
jump ratios of 2:2:4:2:2:2, thus yielding an effec- 
tive resolution of 16384 x 49152 cells. Computations are 
performed using the AMR version of the PLUTO code 
with the HLLC Riemann solver together the spatially 
and tempora lly second-order accu rate MUSCL-Hancock 
scheme. See lMignone et all (|2011[ ) for a detailed descrip- 
tion of the code and implementation methods. 

2.1. Model Parameters and Simulation Cases 

A cylindrical jet equilibrium model is constructed by 
first prescribing radial profiles for density, velocity, mag- 
netic field and then by solving the radial balance momen- 
tum equation for the gas pressure. The details of this 
equilibrium configuration are outlined in the Appendix 
[Al The resulting radial profiles define a family of jet mod- 
els characterized by the hydrogen number density n#, 
longitudinal velocity Vj, temperature T J5 jet to ambient 
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Table 1 

Definition of the simulation sets and 
corresponding parameters. 



Set 


n H 


^ ax (Km/s) 




A(L x ,r, Sv) 


10 4 





152.3 


Ah(L x ,r, Sv) 


10 4 







B{L x ,r, Sv) 


10 4 


10 


422.4 


C(L x ,r,Sv) 


5- 10 4 





152.3 


D(L x ,t, Sv) 


5 • 10 4 


10 


422.4 


E(L x ,r, Sv) 


5 • 10 4 


15 


610.3 



Note. — Different simulation cases are distin- 
guished by the hydrogen density n/j, peak rota- 
tion velocity ^™ ax and magnetic field |£?™ ax | re- 
spectively given in the second, third and fourth 
column. Each set defines a family of models with 
varying X luminosity L x of the central object, 
period and amplitude of the perturbation r and 
Sv. In all simulation cases, the jet radius, tem- 
perature, velocity and density contrast are the 
same and equal to rj = 20 AU, Tj = 2 500K, 
Vj = llOkm/s and 77 = 5, respectively. 

density contrast n = n#/n a and peak rotation velocity 
^max j n present context we restrict our attention 
to purely toroidal configurations and leave models with 
helical magnetic fields (i.e. B z ^ 0) to forthcoming stud- 
ies. Since the ambient temperature is prescribed to be 
T a = 1 000 K, the maximum value of is not a free 
parameter but depends on the rotation velocity. 

Finally, the parameter that controls the degree of pre- 
ionization of the jet material at the base of the jet is the 
X-ray luminosity Lx of the central object, for which the 
ionization at photoionization equilibrium is computed as 
explained in £12.31 

Along with the equilibrium magnetized models we also 
consider purely hydro configurations that, due to an over- 
pressurized beam, cannot establish equilibrium with the 
environment. In this case a conical structure is formed 
during the propagation. 

In the simulations reported here we set the initial jet 
temperature, velocity and density contrast to the values 
Tj = 2500 K, vj = 110 Km/s and 77 = 5, respectively. 
Table [T] summarizes the chosen set of simulation cases 
while we plot in Fig (pQ) the radial profiles for density, 
temperature, velocity and magnetic field. Within each 
set (labeled by a capital letter), the X luminosity of the 
central object, the period and amplitude of the pertur- 
bation are allowed to vary. 

Set A is characterized by no rotations and a relatively 
weak magnetic field and density. As a special case, we 
also include set Ah consisting of purely-HD 2D simula- 
tions. In these cases, the conical expansion favors the 
formation of a decreasing density along the longitudi- 
nal direction. Set B has stronger rotation and (conse- 
quently) magnetic field. Sets C and D are identical to 
A and B (respectively) except that the beam is five time 
heavier. Finally, the last set E has a maximum rota- 
tion velocity ?j^ ax = 15 Km/s and peak magnetic field of 
610 //G. 
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Figure 1. Radial equilibrium profiles for set A (top panel) and 
set B (bottom panel). In each panel we plot density (in 10 4 cm -3 , 
solid line), temperature (in 10 3 K, dotted line), azimuthal velocity 
(in 10 km/s, dashed line) and magnetic field (in 10 -3 Gauss, dash- 
dotted line). 

In previous works on astrophysical jets, we have em- 
ployed a s pecial definition of the in itial perturbation (de- 
scribed in Massaglia et al] (|2005af )). imposing conditions 
that led to the formation of only one shock propagating 
along the jet beam, instead of the usual pair of forward- 
reverse shocks. This approach was preferred because it 
allowed a higher level of control on the energy dissipa- 
tion areas and an easier parallel between the perturba- 
tion parameters and the characteristics of the forming 
Shockwave. 

In the present work however, a time-dependent velocity 
fluctuation is prescribed at the boundary (after a steady 
configuration has been reached) as: 



(4) 



2.2. Initial perturbation 



where r is the period of perturbation (in years). This 
choice is justified by two main reasons: 

(1) the formation of the pair of forward-reverse shocks 
elongates the high intensity line emission area and leads 
to a better agreement with the morphology of the ob- 
served emission knots, and 

(2) our aim of approaching simulation results to obser- 
vational data benefits from less strict conditions on the 
perturbation parameters. 

Moreover, we limit ourselves to three perturbation pe- 
riods since the conditions in which the second and the 
third shock propagate are quite similar. 
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2.3. Pre-ionization fraction 

We analyze the effect of the jet base irradiation by 
X-rays coming from the central TTauri star. Our goal 
is not to model this region in detail, but is limited to 
gain information on reasonable values of the ionization 
of the jet medium at the distance where observations 
and our simulations start, i.e. at r s = 0".l correspond- 
ing to ~ 2 x 10 14 cm. Detailed numerical calculations 
of the combined dynamical, heating-cooling and photo- 
ionization processes in YSO jets are under way and will 
be published in a forthcoming paper. 

Proto-stellar objects show X-ray luminosities 10 28 — 
10 32 ergs s _1 , depending on their mass and pos- 
sibly originating from the magnetized stellar coron a 
(jGlassgold et all l2QQ0b iPreibisch fc Neuhauserl I2005D . 
with possible co ntributions from the jet itself, as dis- 
cussed rece ntly by Skinner 'et al. I (j2011f ) for RY Tau - HH 
938 and bv lGiidel et~al~l (j2011bh for DG Tau. The inter- 
action of a X-ray photon, in the keV energy range, with 
an atom or molecule results in the production of a fast 
photoelectron, the primary, that in tur n generates, col- 
hsionally, a deal of secondary electr ons dGlassgold et al.l 
119970 . We follow the treatment bv lShang et al.l ([2002). 
that ignores the contribution of the primary electrons 
and considers the dominant secondary electrons only. We 
write the energy input Hx by X-rays (energy per unit 
volume per unit time) and the photo-ionization rate £x 
as: 



n (r) f°° 

n * = -f^T J e Lx(E)a pe (E) e- rx y heat dE , 



(5) 



Cx 



Lx(E) 



47rr ■/ £ 6lon 



o- pe (E) e- rx dE . (6) 



In the expression above Lx(E) is the energy dependent 
X-ray luminosity, Eo(= 0-1 keV) is the low-energy cut- 
off, a pe (E) is the cosmic photoelectric absorption cross 
section per H nucleus, yheat is the absorbed fraction of 
the X-ray flux, ei on the energy to make an ion pair, and 
r is the optical pa th in spherical sym metry. Since yheat 
and eion (given bv lShang et al.ll2002[ ) can be considered 
nearly independent of energy, we have 



(7) 



(8) 



Ux = n H (r) y h eat e ion Cx , 
where (jShull fc van Steenberdll985t l 

e ion " 1(H) + I(He) ' 
with 

VH = 0.3908 (1 - ^.4092)1.7592 > 
y He = 0.0554 (1 - ^-4614) 1.666 _ 

In the above relationships 1(H) and I(He) are the ion- 
ization potentials of H and He, x e is the hydrogen ion- 
ization fraction, and 



yheat = 0.9971 [1 - (1 



0.2663\1. 3163i 



r 



specifies the heating fraction. 
The X-ray optical depth tx can be written: 



rx = cj pe (kTx)N, N ■ 



f 

Jo 



n^dr' 



(9) 



where <r pe (E) = a pe (kTx)(keV / E) p and a pe (lkeV) = 
2.27 x 10 -22 cm 2 , kTx = IkeV and the exponent p = 
2.485 is for solar abundances. 

Note that for a thermal spectrum the ionization rate 
(Eq. [6]), becomes 



LxCTpe(kTx) 

47rr 2 e ion 



r 



r p exp[-(£+T X r p ] de. (io) 



where Lx is the total X-ray luminosity and ^ = E/kTx- 
We consider the region close to the inner disk, where 
the disk-wind jet component is originated, and above 
the extended stellar atmosphere, where the stellar-wind 
jet component is be ing launched (see discussion in 
iMatsakos et al.l ((2009)). The medium there is heated and 
ionized by a X-ray flux of luminosity Lx- This region ex- 
tends from a distance r = Rx ~ 10 12 cm (~ 10 R ) from 
the star, i.e. the stellar corona outer radius, up to about 
1 AU, i.e. the inner disk. The radial velocities there are 
small enough and the ionization, recombination, heating 
and cooling timescales fast enough that we can assume 
energetic and ionization/recombination equilibria: 



Ux - £ = o , 
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(Ci + ^)/n-C r (l-/ n )=0, 



(11) 
(12) 



with f n the number fraction of neutral hydrogen atoms, 
^ = ^h(1 — fn + Z) the electron density, nn the total 
hydrogen density and Z (=0.001) the metal abundance 
by number, q, c r are the ionizati on and recombination 
rate c oefficients, respectively (see iDopita fc Sutherland! 
(2003)), and C represents the energy loss term (energy 
per unit volume per unit time) . The loss term is modeled 
accor ding to the SNEq cooling model by iTesileanu et al.l 
(120093) . 

If we assume a very moderate X-ray luminosity Lx — 
10 29 ergs s _1 and a particle density of 10 6 cm -3 , at 
r = Rx we obtain, according to Eqs. [TTIand fT^ equilib- 
rium temperature ~ 12,000 K and ionization fraction 
w 40% close to the jet axis and drops to about 5% 
at 1 AU, at the jet initial lateral border. The ioniza- 
tion/recombination timescales are of the order of months, 
while the heating/cooling ones are about an order of 
magnitude smaller. The matter is then funnelled into 
the jet by dynamical and MHD processes, expands and 
accelerat es reaching velocities of 100 — 200 km s" 1 in a 
few AUs (IZanni et aLll2QQ7t iTzeferacos et al.ll2009f ). One 
may expect a substantial drop in temperature by cool- 
ing, but the ionization fraction, due to long recombina- 
tion timescale, t ~ l/(c r n e ), would remain close to the 
equilibrium one. Thus, the assumption of a residual ion- 
ization fraction in the central spine of the jet of about 
10-20% at 0'M is a quite reasonable one. 

2.4. Post-processing and data analysis 

The output from numerical simulations, that include 
the chemical/ionization network and radiative cooling 
losses, cannot be directly compared with observations. 
Density, velocity and ionization fraction distribution 
must be in fact transformed into surface brightness maps, 
line ratios and Position- Velocity diagrams in a post- 
processing phase. 

The first step in this process is the computation of 2D 
emissivity maps at wavelengths corresponding to atomic 
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transitions of interest, selected by the user. In the 5- 
level atom model considered by the cooling treatment 
implemented in the PLUTO code, there are a few hun- 
dred selectable emission lines. For these computations, 
the ionization state of the matter and the temperature 
in each simulation cell must be known. The simulation 
code PLUTO delivers the detailed ionization state for 
the atomic species H, He, C, N, O, Ne and S. The tem- 
perature is computed from the pressure, density and ion- 
ization state in each cell. 

The second step is the 3D emissivity integration, in 
cylindrical symmetry, done by rotating the 2D emissiv- 
ity maps previously obtained around the z axis. The 3D 
structure is then projected onto a plane perpendicular on 
the line of sight (the emitted power in each emission line 
is integrated over lines parallel to the line of sight), in or- 
der to obtain a surface brightness map similar to the ones 
observed. A simulation of the effects of the PSF of the 
instrument is also added, usually the simulations having 
much higher resolutions than the observational data (in 
order to capture the physics within). The PSF assumes 
a Gaussian form, with user-defined half- width a. For the 
2D surface brightness maps presented in this work, a PSF 
that is roughly 1/4 of the one of HST was employed (HST 
has a resolution of approximately O'M, that means 14 
AU at the distance of Taurus- Aurigae where the sources 
are located). We have chosen to use this smaller PSF 
in order to have, at this stage, a better resolution of the 
output jet structures. In drawing the plots of line ratios 
and surface brightness along the axis of the simulated 
jet, the resolution was reduced to approximately that of 
HST. 

The two steps leading from the PLUTO output data 
to simulated maps of surface brightness are illustrated in 
Fig. [2] 
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100 200 300 400 
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Figure 2. Top-panel: Logarithmic density map from PLUTO 
output; Middle: Emissivity in [SII]6716Ain units of erg cm~ 3 s _1 , 
logarithmic map; Bottom: Surface brightness map in erg 
cm _2 arcsec _2 s _1 , logarithmic, angle jet - LoS 45 deg. 

After the second step of post-processing, a longitudi- 
nal or transversal slit of arbitrary size can be defined on 
the computed surface brightness map, used to compute 
synthetic spectra and position-velocity diagrams. The 
synthetic spectra include the natural and Doppler line 
broadening, and consist of all emission lines selected for 
processing, with customizable spectral range and reso- 
lution. The resulting position-velocity (PV) diagrams 



can be directly compared to the ones derived from ob- 
servations. PV diagrams taken with a slit parallel to 
the jet axis and stepped across the jet or a slit perpen- 
dicular on the jet axis are particularly useful for simu- 
lations that include the rotation of the jet. This is ex- 
pected from models of jet generation, and indications 
of rotation ha ve been detected in several microjets in 
recent works (iBacciotti et all I20021 iWoitas et aLl l2QQ5b 
ICoffev et al.ll2QQ4L l2QQ7h .~ 

It is also possible to extract velocity channel maps in 
custom velocity channels and emission lines, to be com- 
pared with observations. These velocity channel maps 
are of paramount importance in the investigation of jet 
structure. 

3. RESULTS 

We discuss the results of the numerical simulations and 
compare these with observations of emission knots of the 
three sources, for which high-quality observational data 
are available in the literature. The jets obtained with the 
numerical simulations have been projected at an angle of 
45 degrees with the line of sight, taking as a reference 
the case of the RW Aurigae jets. 

3.1. Shocked jet emission 

In the simulation set A, the equilibrium of a cylindrical 
jet is guaranteed by the toroidal component of the mag- 
netic field vector and the density along the jet remains 
uniform, thus the first shock propagates in a constant 
density environment. On the contrary, the second and 
the third shocks in the array travel in the decreasing 
density zone following the propagation of the previous 
shock, ensuring a longer time-span for intense line emis- 
sion. Indeed, following the evolution of the shocks over 
time, one can notice the different behaviour of the second 
shock with respect to the first one, being brighter over a 
larger distance. 

Fig. [3] shows surface brightness maps in three emission 
lines from [SII], [OI], and [Nil] respectively, for a sim- 
ulation type A. In this case the jet variability period 
is 10 years, the perturbation amplitude 50km/s and the 
temperature of the jet material 2 500K (Hydrogen mostly 
neutral before the shock). We can see in this figure the 
sharp decrease of the brightness after the peak of about 
four orders of magnitude over a distance of 50 AU. This 
leaves large dark spaces between emission knots, that are 
not seen in observations. We note that an attempt to al- 
leviate this problem by diminishing the time periodicity 
of the perturbations that evolve in shock waves, lead to 
a decrease in the maximum knot brightness, explained 
by the lower mass flux entering each shock. 

When an X-ray-induced pre-ionization of the pre-shock 
medium is considered (about 19% in Hydrogen), Fig.HJ 
the emission areas behind the shocks are extended com- 
pared to previous case, Fig. [3] (in the figures being pre- 
sented the same moment in the evolution) , and the max- 
imum values of the brightness are higher as well. This 
configuration provides surface brightness maps more sim- 
ilar to observational data, with elongated emission knots 
because of the higher background ratio of ionized ele- 
ments in the jet material. 

Moreover, the presence of pre-ionization leads to the 
increase in the peak surface brightness with factors be- 
tween 2 and 4. This is due to the fact that a pre-existing 
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Figure 3. Simulation in A configuration, no pre-ionization, per- 
turbation amplitude 50km/s, period 10 years, surface brightness in 
[SII]6716A(top), [OI]6300A(middle) and [NII]6583A, units of erg 
cm _2 arcsec _2 s _1 , loglO maps. 
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Figure 4. Simulation in A configuration, pre-ionization 19%, per- 
turbation amplitude 50km/s, period 10 years, surface brightness 
in [SII]6716A(top), [OI]6300A(middle) and [NII]6583A, units of 
erg • cm~ 2 arcsec~ 2 s~ 
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increased number of free electrons fasten the collisional 
ionization and excitation, enhancing the total brightness. 

In Fig. [5] we show a comparison among the simulated 
surface brightness of the shocked jet in [Nil] 6548 Afor 
four different simulation sets (A, B, C, and D from top 



100 

50 r 

T 

-50 ~t 
-100 L 



Surface brightness [Nil] 6548A, PSF 



200 300 
AU 



100 

50 t 

T 

-50 r 
-100 L 



200 300 
AU 



100 I 

so ; 
o ; 
-so ; 

-100 E 



100 



200 300 
AU 



400 



500 



100 
50 t 

r 

-50 ~t 
-100 L 



-14.28 
-15.62 
-1 6.95 
-18.28 
-19.62 
-20.95 
J-22.28 



-14.53 
1-15.87 

-17.20 
-18.53 
-19.87 
1-21.20 
-22.53 



-14.00 
-15.33 
-1 6.67 
1-18.00 
-19.33 
-20.67 
-22.00 



-13.92 
-15.25 
-1 6.58 
-17.92 
-19.25 
-20.58 
J -21 .92 



200 300 
AU 



Figure 5. Surface brightness maps in [NII]6548A, in four sim- 
ulation configurations (A, B, C, and D), with pre-ionization and 



the same set of parameters. Units of erg cm 2 arcsec s 
maps. 



log 



to bottom panels), including the pre-ionization of the 
jet material by X-rays. For consistency, the maps are 
drawn at the same evolutionary stage and the "variable" 
parameters were set to the same values. 

The top panel shows the surface brightness map for 
a simulation in setup A, with rather compact emission 
knots and low-intensity gaps between them. The B sim- 
ulation (second panel from top) includes the jet rota- 
tion with a maximum velocity of 10km • s -1 , and pro- 
duces maximum surface brightness lower than in the 
corresponding A cases, but with a reduced decrease in 
brightness in the regions between two successive emission 
peaks. In case C (third panel in Fig.[5j), the propagation 
of the knots is slightly faster with respect to the pre- 
vious cases, because of the higher density (5 • 10 4 cm -3 
instead of 10 4 cm -3 ). In addition, the maximum value 
of the surface brightness is higher than in the otherwise 
very similar results of case A. The results of case D has 
been obtained setting the jet rotation at 10km ■ s -1 and 
density at 5 • 10 4 cm -3 ), and from Fig.[5j bottom panel, 
we see that the morphology of the line emission is similar 
to the one of case 5, but with higher emission intensities 
due to the increased amount of mass load of the jet. 

The purely hydrodynamic case Ah is characterized by a 
larger lateral expansion, thus both the maximum surface 
brightness and the length of the high-intensity zone result 
lower than in the corresponding MHD cases, so it was 
excluded from the comparison in Fig. [5] The results in 
the E cases were very similar to the ones obtained in the 
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D setup, thus not displayed. 

3.2. Comparison with observations 
3.2.1. Observational constraints 

refer to Hubble Space Telescope (STIS instru- 

t) observations of RW Aurigae jets, [M elniko v etall 

(20(% 



Simulated surface brightness, A 



We 
ment 



for DG Tau, 



Bacciotti et al. 
(2000); and 



(12002! ) and 
for HH30, 



Lavallev-Foua uet et al.l 
Hartigan & Morsel (|2007D . 

In Fig.[6j we show the observed surface brightness 
along the jet axis in the three emission doublets of [OI] 
(6300 Aand 6363 A), [Nil] (6548 Aand 6583 A) and [SII] 
(6716Aand 6731A) for the three sources quoted above. 
Hereafter, where no wavelength is specified, the square 
brackets notation refers to the sum of both lines of the 
respective doublet. 
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Figure 6. Surface brightness in along the jets in units of erg 
cm _3 arcsec _2 s _1 , logarithmic plot on the jet axis from observa- 
tions of RW Aurigae redshifted jet, DG Tau and HH 30. 

One can note the overall higher brightness of the three 
emission doublets for DG Tau, in agreement with both 
the higher Doppler velocities measured for this source 
and the presence of an X-ray e mission discovered by 
the C handra Observatory (e.g. iSchneider fc Schmittl 
(2008)), as possibly indicative stronger shock waves. 
Working in the approximation of optically-thin plasma, 
the higher values for the RW Aurigae redshifted jet with 
respect to HH30 jet, despite the similar flow and shock 
velocities, may be explained by the higher declination 
angle of the former with respect to the line of sight and 
the different toroidal magnetic field strength. 

3.2.2. Surface brightness 

As discussed in the previous section, the surface bright- 
ness variation with distance along the jet differs de- 
pending on the case considered. Without including pre- 
ionization (i.e. with the ionization fraction taken in col- 
lisional equilibrium at 2 500K ahead of the shocks), the 
distribution of surface brightness along the jet (Fig. [7]) 
has variations of many orders of magnitude and lower 
peak values with respect to the pre-ionized cases (and 
much lower than observations). 
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Figure 7. Surface brightness, doublets of [SII], [OI] and [Nil] in 
units of erg cm _3 arcsec _2 s _1 , logarithmic plot on the jet axis from 
simulation type A, with no pre-ionization. 

The rotating jet simulated in configuration D is a good 
candidate for the comparison with observations, the de- 
crease of brightness between the high-intensity being less 
pronounced than in the corresponding non-rotating case 
(A) - see Fig.El 

An important increase in brightness is also impor- 
tant for the comparison with observations - shocks with 
higher- amplitude perturbations (higher than 50 km s _1 ) 
are not likely for "slow" jets such as HH30 and RW Aur, 
so the pre-ionization provides a way of enhancing bright- 
ness without going with the simulations beyond the most 
probable parameter range. 
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Figure 8. Surface brightness, doublets of [SII], [OI] and [Nil] in 
units of erg cm - 3 arcsec — 2 s — 1 , logarithmic plot on the jet axis from 
simulation type D, with pre-ionization. 

The decreasing trend of the peak brightness with the 
traveled distance from the jet origin is visible both in 
simulations and observations: at angular distances larger 
than 2", the decrease is approximately one order of mag- 
nitude (Figs. [6] and [8j). This suggests that the knots 
observed in many jets (e.g. HH 34 and HH 111) at dis- 
tances of a few tens of arcseconds from the source are 
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likely to arise from other mechanisms, i.e. jets shear- 
layer instabilities. 

3.2.3. Line emission ratios 

The line emission ratios are indispensable ingredi- 
ents in methods for deriving the physical parameters of 
space plasmas from observations. In the case of stellar 
jets - the forbidden emission doublets of [SII], [01] and 
[Nil] between 6 and 7QQQAare used ("BE" technique, 
iBacciotti et al.l (|1999l) ) for this purpose. For this reason 
the comparison between the observed and simulated line 
ratios is a powerful method of validation for both the 
numerical code and the correct interpretation of obser- 
vational data. 

In the previous ID analyses we considered the emis- 
sion of a single shock at different times while propagat- 
ing along the jet, instead we are now taking snapshots at 
given times of the whole length of the jet and study the 
behaviour of the line ratios as a function of the longitudi- 
nal coordinate. The high numerical resolution achieved 
thanks to the AMR technique allows us to follow not only 
the values in the emission peaks, but also their evolution 
in the post-shock zone as the gas cools. We draw in Fig. [9] 
the results of the calculations, without pre-ionization, of 
three line ratios of forbidden lines in comparison with the 
observed line ratios (symbols) for the first part of the red- 
shifted jet from the RW Aurigae pair. We see that the 
values of the calculated line ratios approach observations 
only for short distances after the shocks. 
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Figure 9. Line ratios between the three doublets of 
[SII]6716+673lA, [OI] 6300+6363 Aand [NII]6548+6583A, loglO 
scale plot on the jet axis from simulation type A without pre- 
ionization, and observations. 

In Figs. [10] and (Fig.fTT]) we show a simulation from the 
A and B sets, respectively, with pre-ionization included. 
In both cases the behaviour of the calculated line ratios is 
much more consistent with observational data, the vari- 
ations between knots remaining in the observed ranges. 

3.3. Position- Velocity diagrams 

In order to illustrate the distribution in velocities of the 
emitting material, the Position- Velocity (PV) diagrams 
are widely used. A spectrum is generated for each pixel 
along the spectrograph slit, and the results are plotted 
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Figure 10. Line ratios between the three doublets of 
[SII]6716+673lA, [OI] 6300+6363 Aand [NII]6548+6583A, loglO 
scale plot on the jet axis from simulation type A with pre- 
ionization, and observations. 



Observed line ratios, log 1 scale 



= 3 E 




+ 


"■ !, ] 














O 


"S\\ 


5/1 6]/ S 16731 ] 






sim 


Oll/TNII] 






sim 


^sirj/roij 






sim 


= SII6716]/[SII6731] 





1 2 3 4 5 

Distance(arcsec) 

Figure 11. Line ratios between the three doublets of 
[SII] 6716+6731 A, [OI] 6300+6363 Aand [NII]6548+6583A, loglO 
scale plot on the jet axis from simulation type B with pre- 
ionization, and observations. 

in units of surface brightness at a certain wavelength on 
a Position- Velocity map. 

If Fig.[12j the output from the PLUTO post-processing 
routines is shown. The top panel is a surface brightness 
map in one of the lines of the [SII] doublet, with the 
user-defined slit from where the data for the PV-diagram 
will be taken. The bottom panel displays the resulting 
PV diagram, in units of surface brightness. The distri- 
bution of brightness is concentrated to the right half of 
the image, corresponding to positive velocities, due to 
the declination angle between the jet axis and the line of 
sight. The enhanced emission knots can be clearly seen in 
the PV diagram, concentrating around radial velocities 
of 90kms _1 . The inter-knot jet material is distributed 
in a range of velocities between -10 and +70kms _1 . 

The PV diagrams are a powerful tool in the modern 
study of the structure of stellar jets, providing more ac- 
curate information on the velocity distribution of the 
emitting material. By the differences in the radial ve- 
locity and asymmetries between opposite parts of the 
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Figure 12. Surface brightness (top panel) with the defined slit 
;/ .4 wide (at the distance of Taurus), and position- velocity dia- 
gram (bottom) for [SII] 673lA. 



4. CONCLUSIONS AND SUMMARY 

Starting from numerical MHD simulations that include 
ionization network and detailed radiative cooling, we 
have obtained synthetic emission maps of surface bright- 
ness at various wavelengths relevant for observations of 
HH microjets. The comparison with observations was 
not limited to surface brightness (along the jet, inte- 
grated in velocity), we have also tried to match the ob- 
served line ratios for different values of the simulation 
parameters. 

We have shown the crucial role assumed by the pre- 
existing ionization in the jet medium, prior to the pas- 
sage of the shock wave, for the line emission properties of 
the corresponding "knot" . We believe that pre-ionization 
will be a key ingredient in future work. This relatively 
high ionization fraction is likely to come from the X-ray 
photoionization of the atoms at the jet base, being ad- 
vected away with the flow conserving its value because 
of the low the recombination rate. The pre-ionization in- 
creases the number of free electrons in the gas and speeds 
up the processes of ionization and excitation at the pas- 
sage of the shock wave. 

Among the simulations performed during this work, 
the B and D sets, that include a toroidal magnetic field, 
rotation of the jet and pre-ionization, seem to compare 
well with observations. Future analyses will address the 
problem of the contrast between the knots and intra- 
knots brightness, that remains higher than observed, for 
performing simulations aiming to reproduce in greater 
detail the emission features of particular objects, with the 
goal to constrain the jet physical parameters and better 
understand the physical mechanisms at work. Moreover, 
a challenging but potentially insightful investigation will 
be the 3D case, that could address the shock misalign- 
ment. 



j et, their rotation ( predicted by models) can be inferred 
([Coffey et al.l 120071 ). Consequently, as both the spatial 
and spectral resolutions of observational data increased, 
these diagrams were geenrated also from the jet models, 
in order to be compared to the ones derived from ob - 
servations (jCerqueira et al.l [2006: Smith & Rosen 20071). 
Arrays of models were devised ([Kajdic et al.ll2006[ ). 

An interesting study underway, where PV diagrams 
from multiple slits will be employed, focuses on DG-Tau 
and RW Aurigae, in the search for rotation signatures. 
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APPENDIX 

RADIAL EQUILIBRIUM SOLUTION 

The equilibrium solution is constructed by considering the radial force balance between pressure, magnetic and 
centrifugal forces under the assumption v r = B r = 0. The equilibrium condition is expressed through the steady- state 
r-component of the momentum equation, which reads 



dp pv% 1 
dr r 2 



1 d(rB^) 2 dB\ 



dr 



dr 



(Al) 



In the present context, we will ignore the effect of a poloidal field component and simply consider cases with B z = 0. 
Density and longitudinal velocity profiles can be chosen to smoothly match their ambient values for r > Rj while the 
azimuthal component of magnetic field is prescribed by 



Bn 



B,(r) = -^71 



exp [— (r/a) 4 



(A2) 
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where a = 0.9 is the magnetization radius and rj is the jet radius. This choice guarantees that at large radii the 
field becomes essentially force- free whereas close to the axis the electric current J z « — 2B m rj/a 2 is approximately 
constant. A convenient profile for the azimuthal velocity is 




.., /2exp[-(r/a) 4 ] /A N 

v ^ = a J^ l P — » 

where the constant a sets the am ount of rotation and the relative importance of the centrifugal to the Lorentz force. 
With these assumptions Eq. (|A1[) can be integrated giving 

M-H+ ^-'iFf"* . (A4) 

where pj is the jet pressure on the axis. Clearly, when a > B mi the gas pressure increases monotonically with r while 
the opposite is true for a < B m . The condition a = B m yields exact balance between rotations and magnetic forces. 

The actual value of a can be expressed in terms of the maximum rotation velocity v™ ax which, in the limit of 
constant density, becomes 

°«*r(I) 1/4 £v7* (A5) 

Finally, in order to specify the magnetic field strength B m: we note th at, b y assigning the equilibrium ambient 
temperature T a = p a p a ^a/(PakB) (where p a is the ambient density), Eq (|A4[) may be solved for the magnetic field 
strength B m giving 

)2 _ , , 2k B (a\\ (T, T„ 



B ™=°- + 7*k^)*{7r^) (A6) 

where ks is the Boltzmann constant, m a is the atomic mass unit, pj is the jet density, r] = pj/p a is the jet to ambient 
dens ity contrast, pj and p a are the mean molecular weights in the jet and in the ambient medium, respectively. Eq 
(|A6|) immediately shows that, for T a < Tj, the magnetic field has a lower threshold value and its strength always 
increases with rotation. 

REFERENCES 



Bacciotti, F., Eisloffel, J. 1999, A&A, 342, 717 

Bacciotti, F., Ray, T.P., Mundt, R., Eisloffel, J., Solf, J. 2002, ApJ, 576, 222 

Cerqueira, A.H., Velazquez, P.F., Raga, A.C., Vasconcelos, M.J., de Colle, F. 2006, A&A, 448, 231 
Coffey, D., Bacciotti, F., Woitas, J., Ray, T.P., Eisloffel, J. 2004, ApJ, 604, 758 
Coffey, D., Bacciotti, F., Ray, T.P., Eisloffel, J., Woitas, J. 2007, ApJ, 663, 350 
Cox, D., & Raymond, J. 1985, ApJ, 298, 651 

Dopita, M.A., Sutherland, R.S. 2003, "Astrophysics of the diffuse universe", Springer 

Glassgold, A.E., Najita, J., Igea, J. 1997, ApJ, 480, 344 

Giidel, M. 2011a, Bulletin of the American Astronomical Society, Vol. 43 

Giidel, M., Audard, M. Bacciotti, F., & al. 2011b, arXiv: 110 1.2780^ 1 

Glassgold, A.E., Feigelson, E.D., Montmerle, T. 2000, Protostars and Planets IV, Tucson: University of Arizona Press (eds Mannings, V., 

Boss, A.P., Russell, S. S.), p. 429 
Hartigan, P., Hartigan, Raymond, J., Hartmann, L. 1987, ApJ, 316, 323 
Hartigan, P., Morse, J. A., Raymond, J. 1994, ApJ, 436, 125 
Hartigan, P., & Morse, J. 2007, ApJ, 660, 426 
Kajdic, P., Velazquez, P.F., Raga, A.C. 2006, RMxAA, 42, 217 
Lavalley-Fouquet, C, Cabrit, S., Dougados, C. 2000, A&A, 356, L41 
Massaglia, S., Mignone, A., Bodo, G. 2005a, A&A, 442, 549 
Matsakos, T., Massaglia, S., Trussoni, E., & al., A&A, 502, 217 

Melnikov, S.Y., Eisloffel, J., Bacciotti, F., Woitas, J., Ray, T.P. 2009, A&A, 506, 763 

Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C, Ferrari, A. 2007, ApJS, 170, 228 

Mignone, A., Te§ileanu, O., Zanni, C. 2009, "Numerical Modeling of Space Plasma Flows: ASTRONUM-2008" , ASP Conference Series, 
Vol. 406, 105 

Mignone, A., Zanni, C, Tzeferacos, P., et al. 2011, accepted for publication in ApJS 
Preibisch, Th.,& Neuhauser, R. 2005, ApJS, 160, 390 

Raga, A. C, De Colle, F., Kajdic, P., Esquivel, A., Canto, J. 2007, A&A, 465, 879 

Rossi, P., Bodo, G., Massaglia, S., Ferrari, A. 1997, A&A, 321, 672 

Schneider, P.C. & Schmitt, J.H.M.M. 2008, A&A, 488, L13 

Shang, H., Glassgold, A.E., Shu, F.H., Lizano, S. 2002, ApJ, 564, 853 

Shull, J.M., & McKee, C.F. 1979, ApJ, 227, 131 

Shull, J.M., & van Steenberg, M.E. 1985, ApJ, 298, 268 

Skinner, S. L., Audard, M. & Giidel, M. 2011, ApJ, 737, 19 

Smith, M.D., Rosen, A. 2007, MNRAS, 378, 691 

Takami, M., Chrysostomou, A., Ray, T.P., & al. 2004, A&A, 416,213 
Te§ileanu, O., Mignone, A., Massaglia, S. 2008, A&A, 488, 429 

Te§ileanu, O., Massaglia, S., Mignone, A., Bodo, G., Bacciotti, F. 2009a, A&A, 507, 581 



Numerical simulations of radiative proto-stellar jets: pre-ionization from X-rays 

Te§ileanu, O., Mignone, A., Massaglia, S. 2009b, Procs. " Protostellar Jets in Context", T.P. Ray, K. Tsinganos &; M. Stute eds., 

Astrophysics and Space Science Proceedings Series, Springer, p. 447 
Tzeferacos, P., Ferrari, A., Mignone, A., & al. 2009, MNRAS, 400, 820 
Verner, D.A., Yakovlev, D.G. 1995, A&AS, 109, 125 

Woitas, J., Bacciotti, F., Ray, T.P., Marconi, A., Coffey, D., Eisloffel, J. 2005, A&A, 432, 149 
Zanni, C, Ferrari, A., Rosner, R., Bodo, C, Massaglia, S. 2007, A&A, 469, 811 



